# faces.exit.poll.R
#
# Model vote share as a function of faces, competitiveness,
# incumbency, partisanship, and controls.

cat('\nRunning exit poll models ...\n')

library(Hmisc)
library(lmtest)
library(MSBVAR)

load('exit.poll.Rdata')     # Data created in setup.exit.poll.R

########################################################################
#CLUSTERED STANDARD ERROR FUNCTION
##modified from Mahmood Arai, http://people.su.se/~ma/mcluster.R 
########################################################################
cl <- function(dat,fm, cluster) {
    attach(dat, warn.conflicts = F)
    library(sandwich)
    M <- length(unique(cluster))   
    N <- length(cluster)              
    K <- fm$rank                   
    dfc <- (M/(M-1))*((N-1)/(N-K))  
    uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
    vj = uj[is.na(uj[,1])==F,]
    vcovCL <- dfc*sandwich(fm, meat=crossprod(vj)/N)  
    vjs = crossprod(vj)
    coefs = coeftest(fm, vcovCL)
    output = list()
    output$coefficients = coefs
    output$vcov <- vcovCL
    return(output)   
}


############################################################################################
# BASIC EXIT POLL MODELS
############################################################################################

# Save out tex file of regressions
source("lmLatexWriter.R")
listOfLms <- list() # For probits
robustSEs <- list() # For probits
listOfLms2 <- list() # For linear prob models
robustSEs2 <- list() # For linear prob models
model.vars <- c("cook","shares.chal.party","shares.inc.party","face.chal","face.inc",
                "inc.yrsinoffice","I(inc.yrsinoffice^2)","inc.age","I(inc.age^2)")
model.robust <- c("cook","shares.chal.party","shares.inc.party","face.chal","face.inc",
                "inc.yrsinoffice","I(inc.yrsinoffice^2)","inc.age","I(inc.age^2)",
                "chal.expend.log")

#
# House
#

# Basic model
house.formula <- paste("inc.vote ~",paste(model.vars,collapse="+"))            
listOfLms[["House 2004"]] <- glm(as.formula(house.formula),
        family=binomial(link="probit"),data=exit.house)
robustSEs[["House 2004"]] <-  cl(dat=exit.house,
                listOfLms[["House 2004"]],cluster = exit.house$district)$coefficients[,2]
listOfLms2[["House 2004"]] <- lm(as.formula(house.formula),data=exit.house)
robustSEs2[["House 2004"]] <-  cl(dat=exit.house,
                listOfLms2[["House 2004"]],cluster = exit.house$district)$coefficients[,2]
                
# Robust model
house.formula.robust <- paste("inc.vote ~",paste(model.robust,collapse="+"))            
listOfLms[["House 2004 "]] <- glm(as.formula(house.formula.robust),
        family=binomial(link="probit"),data=exit.house)
robustSEs[["House 2004 "]] <-  cl(dat=exit.house,
                listOfLms[["House 2004 "]],cluster = exit.house$district)$coefficients[,2]
listOfLms2[["House 2004 "]] <- lm(as.formula(house.formula.robust),data=exit.house)
robustSEs2[["House 2004 "]] <-  cl(dat=exit.house,
                listOfLms2[["House 2004 "]],cluster = exit.house$district)$coefficients[,2]

#
# Senate
#

# Basic model
senate.formula <- paste("inc.vote ~",paste(model.vars,collapse="+"))
listOfLms[["Senate 1992-2006"]] <- glm(as.formula(senate.formula),
        family=binomial(link="probit"),data=exit.senate)

robustSEs[["Senate 1992-2006"]] <- cl(dat=exit.senate,
            listOfLms[["Senate 1992-2006"]],cluster = exit.senate$state.id)$coefficients[,2]
listOfLms2[["Senate 1992-2006"]] <- lm(as.formula(senate.formula),data=exit.senate)
robustSEs2[["Senate 1992-2006"]] <- cl(dat=exit.senate,
            listOfLms2[["Senate 1992-2006"]],cluster = exit.senate$state.id)$coefficients[,2]

# Robust model
senate.formula.robust <- paste("inc.vote ~",paste(c(model.robust,'statepop','as.factor(year)'),collapse="+"))
listOfLms[["Senate 1992-2006 "]] <- glm(as.formula(senate.formula.robust),
        family=binomial(link="probit"),data=exit.senate)
robustSEs[["Senate 1992-2006 "]] <- cl(dat=exit.senate,
            listOfLms[["Senate 1992-2006 "]],cluster = exit.senate$state.id)$coefficients[,2]
listOfLms2[["Senate 1992-2006 "]] <- lm(as.formula(senate.formula.robust),data=exit.senate)
robustSEs2[["Senate 1992-2006 "]] <- cl(dat=exit.senate,
            listOfLms2[["Senate 1992-2006 "]],cluster = exit.senate$state.id)$coefficients[,2]

##
## Pool House and Senate. Not presented, referenced in footnote.
##
#
## Basic model
#both.formula <- paste("inc.vote ~",paste(c(model.vars,"senate"),collapse="+"))
#listOfLms[["Stacked"]] <- glm(as.formula(both.formula),
#        family=binomial(link="probit"),data=exit.both)
#robustSEs[["Stacked"]] <- cl(dat=exit.both,
#            listOfLms[["Stacked"]],cluster = exit.both$state.id)$coefficients[,2]
#listOfLms2[["Stacked"]] <- lm(as.formula(both.formula),data=exit.both)
#robustSEs2[["Stacked"]] <- cl(dat=exit.both,
#            listOfLms2[["Stacked"]],cluster = exit.both$state.id)$coefficients[,2]
#
## Robust model
#both.formula.robust <- paste("inc.vote ~",paste(c(model.robust,"senate"),collapse="+"))
#listOfLms[["Stacked "]] <- glm(as.formula(both.formula.robust),
#        family=binomial(link="probit"),data=exit.both)
#robustSEs[["Stacked "]] <- cl(dat=exit.both,
#            listOfLms[["Stacked "]],cluster = exit.both$state.id)$coefficients[,2]
#listOfLms2[["Stacked "]] <- lm(as.formula(both.formula.robust),data=exit.both)
#robustSEs2[["Stacked "]] <- cl(dat=exit.both,
#            listOfLms2[["Stacked "]],cluster = exit.both$state.id)$coefficients[,2]

varlabs <- list(shares.chal.party="Respondent Shares Challenger Party",
            shares.inc.party="Respondent Shares Incumbent Party",
            face.chal="Challenger Facial Competence",
            face.inc="Incumbent Facial Competence",
            cook="Cook Incumbent Risk",
            inc.pres.margin="Incumbent Party Pres Vote",
            chal.expend.log="Challenger Expenditures (logged)",
            statepop="State Population (millions)",
            statepop.sq="State Population$^2$",
            inc.yrsinoffice="Incumbent Tenure",
            'I(inc.yrsinoffice^2)'="Tenure Squared",
            inc.age="Incumbent Age",
            'I(inc.age^2)'="Age Squared",
            senate="Senate Fixed Effect")
for (yr in seq(1994,2006,2)) { varlabs[[paste('as.factor(year)',yr,sep='')]] <- paste(yr,'Fixed Effect') }

statsLs=list(aic="AIC")
lmLatexWriter(listOfLms,robustSEs,toFile=F,rounder="%3.3f",
            varlabs=varlabs,statsLs=statsLs)
lmLatexWriter(listOfLms,robustSEs,toFile=T,rounder="%3.3f",
            varlabs=varlabs,statsLs=statsLs,
            outfile="TableAndFiguresOutput/ExitPollProbit.tex")
# Linear probability version
lmLatexWriter(listOfLms2,robustSEs2,toFile=T,rounder="%3.3f",
            varlabs=varlabs,outfile="TableAndFiguresOutput/ExitPollOLS.tex")

if (F) {    # Set to TRUE to write out Stata replication files.
    #
    # Write out files for Stata replication of clustered errors.
    #
    write.csv(exit.house,'stata.house.data.csv',row.names=F)
    write.csv(exit.senate,'stata.senate.data.csv',row.names=F)
    write.csv(exit.both,'stata.both.data.csv',row.names=F)

    zz <- textConnection('foo','w')
    write(paste('cd',gsub('/','\\\\',getwd())),zz)
    write('\n* House',zz)
    write('insheet using stata.house.data.csv, comma clear',zz)
    write(paste("probit incvote",
                paste(sapply(model.vars,gsub,pattern="\\.",replacement=""),
                collapse=" "),
                "if senate == 0, robust cluster(district)\n"),zz)
    write('* Senate',zz)
    write('insheet using stata.senate.data.csv, comma clear',zz)
    write(paste("probit incvote",
                paste(sapply(model.vars,gsub,pattern="\\.",replacement=""),
                collapse=" "),
                "if senate == 1, robust cluster(stateid)\n"),zz)
    write('* Both',zz)
    write('insheet using stata.both.data.csv, comma clear',zz)
    write(paste("probit incvote senate",
                paste(sapply(model.vars,gsub,pattern="\\.",replacement=""),
                collapse=" "),
                ", robust cluster(stateid)\n"),zz)
    write(x=paste(foo,collapse='\n'),file='stata.cluster.do')
    close(zz)
}


############################################################################################
# SIMULATE Inter-quartile range effects
############################################################################################

#
# Simulation parameters
#
set.seed(2009)
num.draws         <- 500
num.sims.per.draw <- 10000
face.quantiles <- c(.25,.75)
results <- list()

##########################
# SENATE
##########################

#
# Set simulation X values
#
face.inc <- median(c(senate.dem$comp,senate.rep$comp),na.rm=T)
inc.yrsinoffice <- mean(exit.senate$inc.yrsinoffice)
inc.age <- mean(exit.senate$inc.age)

xc <- list()    # To store x values for 3 types: challenger copartisans, incumbent
                # copartisans, and independents
xc[["constant"]] = c(1,1,1)
xc[["cook"]] = c(1,1,1)     # Likely for incumbent
xc[["shares.chal.party"]] = c(0,1,0)
xc[["shares.inc.party"]] = c(0,0,1)
xc[["face.chal"]] = rep(0,3)
xc[["face.inc"]] = rep(face.inc,3)
xc[["inc.yrsinoffice"]] = rep(inc.yrsinoffice,3)
xc[["I(inc.yrsinoffice^2)"]] = rep(inc.yrsinoffice^2,3)
xc[["inc.age"]] = rep(inc.age,3)
xc[["I(inc.age^2)"]] = rep(inc.age^2,3)
xc <- data.frame(xc,check.names=F)
xc <- xc[,c('constant',model.vars)]   # Make sure order matches output$coefs
rownames(xc) <- c("Independent","Challenger Copartisan","Incumbent Copartisan")

cat(sprintf(paste("\n==================\n",
        "Estimating effect of moving Senate challenger face from the",
        "\n%sth (%0.2f) to the %sth (%0.2f) percentile...\n",
        "%s draws of %s respondents each.\n",
        "==================\n",sep=''),
    100*face.quantiles[1],
    quantile(c(senate.dem$comp,senate.rep$comp),face.quantiles[1],na.rm=T),
    100*face.quantiles[2],
    quantile(c(senate.dem$comp,senate.rep$comp),face.quantiles[2],na.rm=T),
    prettyNum(num.draws,","),prettyNum(num.sims.per.draw,",") )) ; flush.console()
    
# Model results
output <- cl(dat=exit.senate,fm=listOfLms[["Senate 1992-2006"]],cluster=exit.senate$state.id)

results[['Senate']] <- list()
for (type in rownames(xc)) {    # Loop over voter partisan types
    results[['Senate']][[type]] <- matrix(nrow=num.draws,ncol=2)
    for (i in 1:num.draws) {    # Loop over draws from coefficient distribution
        bsim <- as.vector(rmultnorm(n=1,mu=output$coefficients[,1],
                                    vmat=output$vcov)) # Draw coefficients
        for (quant in 1:length(face.quantiles)){    # Loop over face values to consider
            # Re-set face value
            xc[type,"face.chal"] = quantile(c(senate.dem$comp,senate.rep$comp),
                                        face.quantiles[quant],na.rm=T)
            x <- as.matrix(xc[type,])
            thetasim <- as.numeric(pnorm(x%*%bsim))      # Pr(incumbent vote | xc)
            # Draw num.sims.per.draw from bernoulli given thetasim probability
            ysim <- ifelse(runif(num.sims.per.draw)<thetasim,1,0)
            results[['Senate']][[type]][i,quant]<- mean(ysim)       # Store expected values
        }
    }
    results[['Senate']][[type]] <- as.data.frame(results[['Senate']][[type]])
    colnames(results[['Senate']][[type]]) <- paste('quantile',face.quantiles,sep='')
    results[['Senate']][[type]]$effect <- results[['Senate']][[type]][,2]-results[['Senate']][[type]][,1]

    cat( sprintf("%s voter effect is: %1.3f [%1.3f,%1.3f].\n",
            type,mean(results[['Senate']][[type]]$effect),quantile(results[['Senate']][[type]]$effect,c(.025)),
            quantile(results[['Senate']][[type]]$effect,c(.975))) )
}

##########################
# HOUSE
##########################

#
# Set simulation X values
#
face.inc <- with(aeh.house,median(ifelse(dem.incum==1,dem_comp,
                                ifelse(rep.incum==1,rep_comp,NA)),na.rm=T))
inc.yrsinoffice <- mean(2004-aeh.house$IncFirstYrInOffice)
inc.age <- mean(2004-aeh.house$IncBirthYr)

xc <- list()    # To store x values for 3 types: challenger copartisans, incumbent
                # copartisans, and independents
xc[["constant"]] = c(1,1,1)
xc[["cook"]] = c(1,1,1)     # Likely for incumbent
xc[["shares.chal.party"]] = c(0,1,0)
xc[["shares.inc.party"]] = c(0,0,1)
xc[["face.chal"]] = rep(0,3)
xc[["face.inc"]] = rep(face.inc,3)
xc[["inc.yrsinoffice"]] = rep(inc.yrsinoffice,3)
xc[["I(inc.yrsinoffice^2)"]] = rep(inc.yrsinoffice^2,3)
xc[["inc.age"]] = rep(inc.age,3)
xc[["I(inc.age^2)"]] = rep(inc.age^2,3)
xc <- data.frame(xc,check.names=F)
xc <- xc[,c('constant',model.vars)]   # Make sure order matches output$coefs
rownames(xc) <- c("Independent","Challenger Copartisan","Incumbent Copartisan")

cat(sprintf(paste("\n==================\n",
        "Estimating effect of moving House challenger face from the",
        "\n%sth (%0.2f) to the %sth (%0.2f) percentile...\n",
        "%s draws of %s respondents each.\n",
        "==================\n",sep=''),
    100*face.quantiles[1],
    quantile(c(aeh.house$dem_comp,aeh.house$rep_comp),face.quantiles[1],na.rm=T),
    100*face.quantiles[2],
    quantile(c(aeh.house$dem_comp,aeh.house$rep_comp),face.quantiles[2],na.rm=T),
    prettyNum(num.draws,","),prettyNum(num.sims.per.draw,",") )) ; flush.console()
    
# Model results
output <- cl(dat=exit.house,listOfLms[["House 2004"]],cluster = exit.house$district)

results[['House']] <- list()
for (type in rownames(xc)) {    # Loop over voter partisan types
    results[['House']][[type]] <- matrix(nrow=num.draws,ncol=2)
    for (i in 1:num.draws) {    # Loop over draws from coefficient distribution
        bsim <- as.vector(rmultnorm(n=1,mu=output$coefficients[,1],
                                    vmat=output$vcov)) # Draw coefficients
        for (quant in 1:length(face.quantiles)){    # Loop over face values to consider
            # Re-set face value
            xc[type,"face.chal"] = quantile(c(aeh.house$dem_comp,aeh.house$rep_comp),
                                        face.quantiles[quant],na.rm=T)
            x <- as.matrix(xc[type,])
            thetasim <- as.numeric(pnorm(x%*%bsim))      # Pr(incumbent vote | xc)
            # Draw num.sims.per.draw from bernoulli given thetasim probability
            ysim <- ifelse(runif(num.sims.per.draw)<thetasim,1,0)
            results[['House']][[type]][i,quant]<- mean(ysim)       # Store expected values
        }
    }
    results[['House']][[type]] <- as.data.frame(results[['House']][[type]])
    colnames(results[['House']][[type]]) <- paste('quantile',face.quantiles,sep='')
    results[['House']][[type]]$effect <- results[['House']][[type]][,2]-results[['House']][[type]][,1]

    cat( sprintf("%s voter effect is: %1.3f [%1.3f,%1.3f].\n",
            type,mean(results[['House']][[type]]$effect),quantile(results[['House']][[type]]$effect,c(.025)),
            quantile(results[['House']][[type]]$effect,c(.975))) )
}
cat('\n\n')

#
# First differences plot.
#
if (toEPS) postscript('TableAndFiguresOutput/FaceIQREffects.eps',
                        width=6,height=6,horizontal=F)
par(mar=c(5.1,4.1,1.1,1.1))
ylim <- c(-0.1,0.1)
plot(x=c(1,(1+2*nrow(xc))),y=ylim,type='n',ann=F,axes=F)
axis(1,at=c(1:nrow(xc),(nrow(xc)+2):(2*nrow(xc)+1)),
        rep(gsub(' ','\n',rownames(xc)),2),cex.axis=.8,las=2)
axis(2,at=seq(from=ylim[1],to=ylim[2],length.out=5),
    round(100*seq(from=ylim[1],to=ylim[2],length.out=5),2),las=2)
title(ylab='Percentage Point Effect on Challenger Vote')
xloc <- 1
abline(h=0,lty=2,col='gray')
for (chamber in c('Senate','House')) {
    text(x=xloc,y=ylim[2],paste(chamber,'Effects'),pos=4,cex=1)
    for (type in rownames(xc)) {
        # Flip direction of effects from pro-incumbent to pro-challenger
        res.vect <- -1*results[[chamber]][[type]]$effect
        points(x=xloc,y=mean(res.vect),pch=19,cex=2)
        lines(x=rep(xloc,2),y=c(quantile(res.vect,.025),quantile(res.vect,.975)),
                lwd=2 )
        xloc <- xloc + 1
    }
    xloc <- xloc + 1
}
if (toEPS) dev.off()
